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Systems with purely off-diagonal disorder have peculiar features such as the localization- 
delocalization transition and long-range correlations in their wavefunctions. To motivate possible 
experimental studies of the physics of off-diagonal disorder, we study in detail disordered discrete¬ 
time quantum walk in a finite chain, where the diagonal disorder can be set to zero by construction. 
Starting from a transfer matrix approach, we show, both theoretically and computationally, that 
the dynamics of the quantum walk with disorder manifests all the main features of off-diagonal 
disorder. We also propose how to prepare a remarkable delocalized zero-mode from a localized and 
easy-to-prepare initial state using an adiabatic protocol that increases the disorder strength slowly. 
Numerical experiments are also performed with encouraging results. 
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I. INTRODUCTION 

Quantum walk (QW) has been a subject of great the¬ 
oretical and experimental interests. Among many QW 
protocols, discrete-time QW is the simplestwhere it can 
be seen clearly how QW can differ strongly from classical 
random walk due to quantum interference effects. For 
example, an initially localized state in QW will spread 
ballistically, which is much faster than classical random 
walk whose mean square displacement is proportional to 
time. Due to this feature, one potential application of 
QW models is towards a fast search algorithm ^ in quan¬ 
tum computation'^. As a very recent direction, QW is 
shown to be useful in understanding topological phases 
of matter in periodically driven systems*’"’. 

On the experimental side, two early QW experiments 
in 2005 used either linear optical elements® or nuclear- 
magnetic resonance systems^. Since 2007, a variety 
of physical systems has been exploited to realize QW, 
including trapped ions®’®, trapped atoms in a spin- 
dependent optical lattice*®, photons in an optical waveg¬ 
uide array’and photonic walks with interferome¬ 
ters*®^**'. Very recently, a photonic quantum walk with¬ 
out interferometers was realized*®, in which photons walk 
in the orbital angular momentum space. 

The topic of this work is on QW in the presence of 
some disorder. Previously, it was numerically found that 
some behavior of disordered QW seems to reflect the 
physics of off-diagonal disorder (ODD)*® in condensed- 
matter physics. The so-called ODD was first noticed 
in studies of one-dimensional (ID) tight-binding mod¬ 
els (TBMs) with random hopping potential and constant 
on-site potentiaP®’^*. Compared with the more famil¬ 
iar disorder model where the on-site potential (diagonal 
term in the lattice-site representation) is random but the 
hopping is constant, ODD leads to peculiar physics, such 
as delocalization at zero energy, power-law wavefunction 
correlation, and so on^®^^®. Specifically, the localization 
length £{u}) in ID TBM with pure ODD is related to 


energy oj via 

£{u}) (X I Intcj. (1) 

As the energy oj approaches 0, the localization length £ 
diverges, indicating a delocalization transition at w = 0. 
At the same time, singularity in the density of states 
(DOS) emerges at w = 0, with the explicit DOS expres¬ 
sion given by 

p{uj) oc |wln®a;|“*. (2) 

Furthermore, the delocalized eigenstate has an unusual 
long-range correlation. It is shown that its ensemble aver¬ 
aged two-point correlation decays polynomially with the 
exponent —3/2 under the condition of strong disorder 
and large two-point separation’^®’’^*’®®. It was pointed out 
earlier that this is a manifestation of the actual stretched 
exponential-decay profile of the wave function'® ’“®‘*, i.e., 
ip{x) oc exp(—7|a: — ccol*/*^), where 7 is a constant. One 
may naively say that a wavefunction like this is quite lo¬ 
calized. However, its Lyapunov exponent is apparently 
zero (which indicates that the state is delocalized®*) be¬ 
cause there is no exponential localization behavior. 

As we have learnt from decades of studies, quite a few 
theoretical models with disorder can be used to man¬ 
ifest and digest the physics of ODD. Such models in¬ 
clude a special disordered linear chain of harmonic oscil¬ 
lators investigated by Dyson®^’®®’®®, a ID Dirac model 
with random mass and some types of disordered ID 
spin chains^®’^*’®®, 2D Dirac fermions subject to a ran¬ 
dom vector potential", a ID random hopping model 
consisting of several parallel bipartite sublattices'*®, sys¬ 
tems with correlated off-diagonal disorder ’®’®® or random 
long-range hopping®*, and graphene with ODD*®. In 
contrast to these theoretical developments, experimen¬ 
tal progresses on the physics of ODD have been rather 
limited. Doped CuGeOs is effectively a disordered spin- 
Peierls system possessing ODD*®“®®. There phenomena 
like phase transitions and long-range orderings were be¬ 
lieved to be related to the physics of ODD. However, di¬ 
rect observation of physical properties like the correlation 
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exponent —3/2 was not possible in such a system. Other 
than spin-chain realizations, few experiments concerning 
ODD were reported. We note a possible experimental 
approach based on cold atoms under the so-called tri¬ 
pod scheme^®’^®, but the actual experiment has not been 
done. Only very recently, Keil et al demonstrated that 
a chain of optical waveguides could be used to realize an 
effective ID Dirac model with random mass^^. In partic¬ 
ular, with coupled series of optical chains, the authors of 
Ref.'^^ observed the long range correlation (in a certain 
range) characterized by the correlation exponent —3/2. 

To motivate more possible experimental studies of 
ODD models and to demonstrate one more promising 
application of QW, we consider in this work a discrete¬ 
time QW in a finite chain (for simplicity we refer to it 
as “QW” throughout the paper) and reveal theoretically 
how this problem is closely connected with the issue of 
ODD. Our work is inspired by an early numerical study 
by Obuse and Kawakami'*^, which showed clear signa¬ 
tures of the physics of ODD in disordered QW. Specifi¬ 
cally, we first analytically demonstrate the explicit con¬ 
nection between a TBM with ODD and disordered QW. 
In so doing we focus on a specific delocalization transition 
energy, the zero quasi-energy, which was also considered 
in Ref.We then show how some simple adiabatic pro¬ 
tocols, starting from an exponentially localized 0-mode 
(i.e., the 0 quasi-energy eigenstate), can be converted to 
a peculiar 0-mode possessing the physics of ODD, with 
satisfactory fidelity and relatively short duration of the 
protocol. As such, we may make use of some existing 
QW experimental set-ups to observe the unique physics 
of ODD. Indeed, our numerical experiments indicate that 
the results agree with theoretical predictions very well, 
including the —3/2 correlation exponent. One advantage 
of this QW approach is that the diagonal disorder does 
not exist by construction, so that the results are free of 
any possible contamination due to diagonal disorder. 

This paper is organized as follows. In Sec. II, we will 
introduce a model of disordered QW in a finite chain. 
Analysis of the model is based on the transfer matrix 
formalism. Sec. Ill is devoted to some formal connec¬ 
tions between our QW model and a TBM with ODD. In 
Sec. IV we shall focus on the preparation of special states 
that best manifest the peculiarities of ODD. The associ¬ 
ated results from our numerical experiments will be also 
presented and discussed. Sec. V concludes this work. 


II. DISORDERED QW IN A FINITE CHAIN 


The standard discrete-time QW is defined via a single 
particle with two internal degrees of freedom. For con¬ 
venience, we refer to its internal states as “spin-up” and 
“spin-down”. The QW protocol consists of two opera¬ 
tions, a rotation of spin through operator R, followed by 
a shift operation by S. Without loss of generality, we 
consider a rotation around y axis by an angle 20, such 


that R = e : 


m = 


f cos 9 
\^sin0 


— sin0\ 
cos 6 J ■ 


( 3 ) 


The operator R rotates the spin at each site, and then the 
spin-up component walks to the right, whereas the spin- 
down component walks the left. Such spin-dependent 
shift operation is implemented via the operator S\ 


OO 

5'= X! (!”- +1) (”■! ® It) (tl +1) (»T-I14^) (4^1) ■ 

n— — oo 

( 4 ) 

The overall one-step quantum walk operator (without 
disorder) is then given by 


Ddt = S 



( 5 ) 


The above described QW can be restricted to a finite 
regime^®’®'*’®'^ through total-reflection coin operators R± 
at two boundaries, with R± defined as 


R± 




/cos(±f) -sin(±f)\ 
Vsin(±f) cos(±f) ) ■ 


Note that R± preserves the particle-hole symmetry and 
conserves the probability inside a finite QW chain. R± 
turns spin-down to spin-up, and vice versa. Since the coin 
operators at two boundaries can be either i?_|_ or R _, we 
could have 4 choices of boundaries as [i?(0o)) ^(6*Af-i-i)] = 
{R±,R±). In the following we mainly choose (i?_,i?+) 
as our boundary condition. Studies of other boundary 
conditions can be found in Appendix B. As depicted in 
Fig. 1, our QW model has totally N + 2 sites, with N of 
them being bulk sites. 

Next we introduce disorder to the QW model, by con¬ 
sidering a perturbation to the local rotation angles 0„, 
i.e.. 


0„ = 0-H(5„ forn= 1, 2, ... TV. (7) 

Here 6 is identical for different sites n, while Sn G [~A, A] 
may differ from site to site, giving rise to a disordered 
QW on a finite number of sites. 

For such a finite-site QW system with a disordered 
bulk specified by 0„, we can still define a mapping oper¬ 
ator U, which can be adpated from the C/dt in Eq. (5) 
[that is, R{9) representation of dif¬ 

ferent QW sites, U can be expressed explicitly as a 
2{N -I- 2) X 2{N + 2) matrix. As a mapping operator, 
U is unitary with eigenvalue 6*“^: 

C/|V-) = e*^|V'), (8) 

where uj is the quasi-energy eigenvalue of U, \ip) is the 
associated eigenstate characterized by 

IV') = (Q^O Po Oil ■ • ■ OiM+l Pn+i)^, 


(9) 
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FIG. 1. (Color online) Set-up of our finite-chain QW with 
disorder, with totally -|- 2 sites, where site 0 and -|- 1 
are the boundary sites with reflection operators R- and 
Rotation operators of bulk sites with n = 1,2 ... N — 1, N 
depend on the local angle 9n, which fluctuate from site to site. 
The red slashes connect spin components /3„ and On+i, as 
they form the new “spinor” in our transfer matrix formalism 
elaborated in our main text. 

with (■••)"'" being the transpose operation. Because of 
the special choices of rotation operators at two bound¬ 
aries, the first and last rows, and the first and last 
columns of U have entries 0 only. Upon removing these 
rows and columns, U becomes a 2{N -|- 1) x 2{N + 1) 
matrix. Correspondingly, the entries ao and /3n+i in the 
eigenstate \iIj) can be also removed. 


which further reduce to 

(^A=r ( ^ /^sin 6 »jv-si\ 

°(^-sin6»oj’ Vajv+ij ^ 6*“" )’ 

(14) 

Using the boundary conditions in Eq. (14) , the chain 
relation in Eq. (11), as well as Oq = —tt/2 and 9n+i = 
'n 12, we finally obtain the following equation that carries 
all the information of Eq. (8): 

CN = Tj, ■ Tm-1 • • • • T2 • Ti • Co . (15) 

For a specific realization of disorder, only particular val¬ 
ues of the quasi-energy oj satisfy Eq. (15). The coeffi¬ 
cients cn and Cg can be determined from Eq. (15) and 
the normalization of lip). 

To conclude, the TMF reduces a matrix equation with 
dimension 2(A^-1-1) x 2(Af-|-l) [Eq. (8)] to a chained matrix 
equation connecting N matrices, each of dimension 2x2 
[Eq. (15)]. This framework will be used later. Indeed, in 
the following we will not return to the original Eq. (8) 
but just focus on Eq. (15). 


A. Transfer matrix formalism 

In solving Eq. (8), one obtains the following recursive 
relation between the entries of the eigenstate lip). 

J O^jiC — CXn—i cos 0 ^ 1 — 1 Pn— 1 sfo On—i, ( 10 ) 
= a„+isin6»„+i-b/3„+iCOS0„+i, 

with n S [1,./V]. Such relations can be expressed in the 
following matrix form: 


Pn \ _ rp (f^n—1 

Q^n+1 


with 


T„. = 


sec On — tan On 
— tan On 6“*“^ sec Or 


( 11 ) 


( 12 ) 


Here is the transfer matrix^*^ at site n. In Eq. (11), 
the neighboring spinors’ components /3n-i and form 
the new “spinors” (See Fig. 1), and they are chained 
through local transfer matrices. Disordered parameter 
On and quasi-energy uj are contained in these matrices. 
This allows us to deal with disorder explicitly. This is 
one known advantage of the transfer matrix formalism 
(TMF)-''’2.53. 

Given the chain relation between entries of the eigen¬ 
state \'ip) in Eq. (11), we still need to handle the boundary 

situations with care, i.e., 

n in Eq. (10) to be 0 and N, we obtain 



: ao cos 00 — I5q sin Oq, 

■ cx-N+i sin 0Ar+i -|- Pn+i cos 0Ar+i, 


(13) 


B. Special quasi-energies and the implication of 
ODD 


By observing the transfer matrix in Eq. (12), we no¬ 
tice that UJ = 0,±7r/2,7r are special quasi-energies. For 
example, when a; = 0, the transfer matrix reduces to 


Tn = sec On -1 — tan t 


(16) 


where I is the identity 2x2 matrix. Such simple transfer 
matrices can be exactly diagonalized in the basis of ctj, , so 
that the product of all the transfer matrices can be eas¬ 
ily calculated. This being the case, whether uj = 0,±7r/2, 
or TT satisfies Eq. (15) can be checked without difficulty. 
If UJ is not equal to one of these special values, then it 
is virtually impossible to analytically check Eq. (15) be¬ 
cause the product of these transfer matrices is hard to 
evaluate. 

If UJ assumes one of these special values, the corre¬ 
sponding eigenstates can be also analyzed in a straight¬ 
forward manner. Take again the case of w = 0 as an 
example. When w = 0, from Eq. (16) we get 


N 


— 2 ('^+ + ^-)I + 2 (■^+ ~ 

n—1 


N 


(17) 


with A+ = A_^ = tan ( ^ ^ . 


And the “spinors” at both ends of ['(/') are proportional to 
1^ 


, i.e., the eigenvector of Ux, obtained from Eq. (14). 


Substituting Eq. (17) into Eq. (15), we get 


Co , /I 
— - 

Cn 


(18) 
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which obviously holds by an appropriate choice of cq/cn- 
Therefore, w = 0 is indeed a quasi-energy solution of the 
disordered QW system. 

In Eq. (17), if fluctuates around 0 or tt (i.e., 0 = 0 or 
tt) , In IA+1 will follow unbiased diffusion process around 
0, so |A+| « 1 for large N, which means that exponential 
decay of the eigenstate |'0) does not occur. This quanti¬ 
tative analysis resembles that of off-diagonal disordered 
TBM'-^'’’^^ so we suspect that our model also displays 
the physics of ODD. Indeed, later in Sec. Ill we shall 
show that a; = 0 is the localization-delocalization transi¬ 
tion quasi-energy, and Dyson’s singularity emerges there, 
provided that 9n takes values randomly from a box dis¬ 
tribution [—A, A]. If On fluctuates around values other 
than 0 or TT, |A+| will increase or decrease exponentially, 
resulting in the localized 0- or 7r-mode, which we believe, 
is related to those topologically protected edge states cur¬ 
rently being studied^*^. 

In the rest of this paper, we focus on the quasi-energy 
a; = 0 and quasi-energies in its vicinity. In Appendix C, 
we shall discuss those cases with quasi-energy values 
other than 0 or tt. 


III. PHYSICS OF ODD 

As introduced in Sec. I, ODD is quite different from 
diagonal disorder and leads to peculiar properties. For 
our QW model, here we attempt to derive its DOS and 
localization length, keeping mind that it is possible for a 
delocalization transition to occur at some special quasi¬ 
energy values. 


A. Analyzing quasi-energy values 


We start with Eq. (15) by considering its alternative 
form after some transformations: 



/ COSLO 

i sina; \ 


(A 

cl . . 

V 2 sincj 

cos a; ) 

p ■ 

(o) 


with 


N 


^=n 

n—l 


/ tan dn 

0 \ 

/ COS a; 

i sina;\ 

[1 0 

cot 

\i sin a; 

cos a; J 


(19) 


where = f — • The detailed derivation can be found 

in Appendix A. Note that if and only if w takes the actual 
quasi-energy value, then Eq. (19) will be satisfied. In 
particular, it is now obvious to observe from Eq. (19) 
that w = 0 is one quasi-energy value. To derive DOS, 
we need to analyze other quasi-energy values allowed by 
Eq. (19). To that end we first re-interpret Eq. (19), which 
is inspired by Schmidt’s work ' that treats spinors linked 
by transfer matrices as vectors in a plane. 

Let us consider a complex plane with x-axis denoting 
the real part, while y-axis denoting the imaginary part. 

In Eq. (19), the initial “spinor” ( 3, ) can be treated as 


a vector lying in the real axis with length 1 pointing in 
the positive direction. So from now on, we refer to the 
“spinor” as a “vector”. Let 


R = 


cos uj I sm oj 
i sin w cos uj 


and Cn = 


tan dn 0 
0 cot -dn 


(20) 

so R and Cn do the job of P in Eq. (19). Consider a 

vector Vn = ( ). Its angle with respect to positive 

\*2/n / 

cc-axis is (/)„, and tan^„ = yn/xn- According to Eq. (19), 
we define 


^n+l — Cn * R ■ V^, 


( 21 ) 


with n = 1, 2, • • • N, and Vi = . Hence, we can in¬ 

terpret Eq. (21) (and Eq. (19) thereafter) as the following 
(see also Fig. 2): R rotates vector v„ counter-clockwise 
by an angle w, followed by stretching in cc-coordinate by a 
factor tan and y-coordinate by the factor cot "On (due 
to Cn), and then v„_|_i is reached with the following re¬ 
lation 


tan (/)„+! = tan (0„ -|- uj) cot^ (22) 

In Eq. (19), the initial vector and final vector v/ 

are both j 'V'i = vi, Vf = R ■ vjy+i, so tan 01 = 

tan(0Ar+i -I- w) =0. As such, Eq. (19) presents such a 
physical picture: a vector initially located in positive x- 
axis is rotated and stretched or contracted, repeatedly, 
and after a final rotation, it lands back on the x-axis. 
Therefore, 


01V- 1-1 +u} = jn. (23) 

Note that uo has the period of 27r, so we assume w S 
[—7r,7r]. Through interpreting Eq. (19) this way, we are 
now ready to derive the DOS near uj = 0. Without loss 
of generality, we consider a small positive quasi-energy 

UJ. 

Regarding the rotating and stretching and contracting 
processes, there are two important factors to be noted. 
First, 0„ does not increase monotonically with respect 
to n. 0„+i could be smaller than 0„ (see Fig. 2). How¬ 
ever, 0„ has a tendency to increase because the positive uj 
forces v„ to rotate counterclock-wise. Besides, a vector 
v„ can never cross x and y-axis clockwise. For exam¬ 
ple, if v„ is inside the first quadrant, then tan0ji and 
cot^ are positive, so for tan0„_|_i in Eq. (22) to be 
negative (i.e., crossing the axis), tan(0„ -|- uj) must be 
negative. Therefore, only the rotation R can bring a vec¬ 
tor from one quadrant to another, while the stretching 
and contracting operation C„ cannot. The vector v„ can 
only drift away by crossing the positive y-axis. Thus, in 
Eq. (23), j is always a positive integer. Second, in a single 
realization of disorder, the following equation holds 

01V+l(Wf,) > (pJV+l(uJa) for UJb > UJa- (24) 
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monotonous function of w. This feature is checked in our 
numerical studies. 

Given the two factors above, we can now count the 
number of states between quasi-energies 0 and ui. Sup¬ 
pose that the corresponding vector of a; sweeps an angle 
in-between jtt and (j -I- l)7r, then there exists j quasi¬ 
energies wi • • • Wj that are the solution of the systems, 
and their vectors sweep angles tt - • • jtt correspondingly. 
Therefore, the number of states between 0 and uj is j, 
and specifically, 


If j < 


0Af(w) 


• U) 


< j +1, 


(26) 


and 


UJi < UJ2 < ■ ■ ■ < OJk ■ ■ ■ < ^j-l < < w, 

with k = 

TT 


(27) 


Here k G [1, j] and it is an integer. Next, we derive the 
integrated DOS from the total number of states. 


B. Integrated density of states 


The general form of the integrated DOS normalized 
over the number of sites is 

/ to 

p{uj')duj'. (28) 

-OO 


FIG. 2. (Color online) The operations in Eq. (19) illustrated 
via a complex plane with thex-axis denoting the real part 
(the first component of the spinor) and the y-axis denoting 
the imaginary part (the second component of the spinor). In 
the first quadrant, from top to bottom, the four vectors are 
v„+i, Vmid, vjj+iand v„. C„ ■ R acts on v„ to jet v„+i (if 
contracted) or (if stretched). Specifically, R rotates v„ 
by angle oj to get Vmid; then C„ will stretch or contract v^id, 
In panel (a), Vmid’s angle is less than 7r/4, while in panel (6) 
its angle is larger than 7r/4. Hence the length of v„+i in panel 
(a) is smaller than in panel (6), whereas the opposite is true 
for V„+i. 

To prove this relation, we show that given <()„ > (j)'^ and 
UJ > uj', then <j)n+i > 4>n+i- We assume that (pn and (f>'^ 
are quite close and and within the same quadrant, say 
the first quadrant. Then it is easy to see that 

tan^„_i_i — tan0„_|_]^ 

= [tan(())„ -I- w) — ta.n{(j)'^ -|- uj')] tan^ > 0, 

so we get (pn+i > (Pn+i- This conclusion can be easily 
proved in other quadrants, too. Hence, starting with 
the same initial condition pi = 0 and same realization 
of disorder, after N cycles, the associated (j)N+i{uj) is a 


Here p{uj) is the density of state (DOS). In QW, particle- 
hole symmetry is present^*’, so quasi-energy uj is symmet¬ 
ric with respect to 0. There are an equal number of posi¬ 
tive and negative quasi-energy states so that Nj{0) = 0.5. 

As shown in the previous section, the total number of 
states between quasi-energies 0 and ui is j, and 

j = + w)/7r], (29) 

where [x] denotes the largest integer less or equal to x. 
So in our case, 

iV,(a;)-iVz(0) = ^^. (30) 

Now we need to evaluate j. 

As shown in Eq. (21), v„+i can be obtained from v„ 
after the operation Cn ■ R. The initial vector will ex¬ 
perience totally A^-|-l operations to reach the final vector 
"Vf. To see this, we add a matrix Cn+i with dAr+i = 0 
to the right of Eq. (19). It is the identity matrix so that 
Eq. (19) holds. From to v/, the vector has passed 
many quadrants. We can define Nq to be the number of 
operations required for the vector to leave the q-th quad¬ 
rant since entering it. Obviously, the summation of all 
the Nq equals to Af -f 1: ^ Afq = Af -|- 1. 
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From Vi to v/, the vector rotates totally by an angle 
about jn after N + 1 operations (see Eq. (29)) so the 
number of quadrants passed is 2j and 

/ 1 \ _ 

= = ^ + ( 31 ) 

q=l 9=1 / 

Hence, we have this formula^*, 

iV,(a;)-7Vr(0)=j7(iV + l) = ^, (32) 

and Ng is the average number of operations required to 
pass one quadrant since entering it. Equation (32) resem¬ 
bles Eq. (21) in the paper by Eggarter and Riedinger^^. 
Though we approach the DOS through counting the 
number of states like what was done in Ref.^^ we are 
able to achieve this step by first introducing the trans¬ 
fer matrix approach when analyzing the spinors in our 
QW model. More importantly, because the above ex¬ 
pression for counting the number of states is similar to 
that in Ref.^', we can now analogously derive the DOS 
near a; = 0. 


C. Derivation of the DOS 

In the previous subsection, the integrated DOS is de¬ 
rived in Eq. (32), but with one parameter Ng to be de¬ 
termined (which represents the average number of oper¬ 
ations required to pass one quadrant). Without loss of 


generality, we consider the first quadrant. 

Let Zn = cot (j)n- From Eq. (22) we have 


I - (tana;)/ 2 ;„ 2 q 

-2^71+1 — , tan Un- 

1 Zn tana; 

(33) 

We define = In for Zn 0 or oo. When 


tanw Zn (tana;)”^. 

(34) 

one approximately has 


Un+I « u„ -I- In (tan^ 'd„) . 

(35) 


Since '&n is taken randomly from this interval [7r/4 — 
A,7r/4 -I- A], we can conclude that Un executes a ran¬ 
dom walk^i. One may notice that the fraction factor in 
Eq. (33) is always smaller than 1 for positive so the 
random walk in Eq. (35) is accompanied with a small 
negative drift. However, if the vector falls in the second 
quadrant, the fraction factor will be always larger than 
1, such that the random walk has a small positive drift. 
The two drifts cancel each other approximately. 

When Un approaches the endpoints of the interval 
in (34), the approximation in (35) no longer holds. Here 
we analyze the situations upon approaching the end¬ 
points to show that they are similar to the situations 
analyzed in Ref.^ft If this is true, then the derivation 
there can be adopted here without much modification. 


For Zn ~ (tana;)“^ (approaching the large Zn limit), 
then Zn+i « (l/2)z„ tan^ according to Eq. (33). The 
net shrinking factor (1/2) in this expression indicates 
that Zn+i will not keep growing. So Umax = — In tan w 
can be considered as the reflection barrier as in Ref.^^. 
We can also view the reflection as the manifestation 
that the vector can never cross x-axis clockwise (see 
Sec. HIA.). 

In the other extreme where z„ « tanw (approach¬ 
ing the small limit), the numerator in Eq. (33) will 
be much smaller than 1 so that Zn+i << Zn, indicat¬ 
ing a sharp decrease in z„. Once gets slightly be¬ 
low tanw, Zn+i will be negative, indicating that the vec¬ 
tor moves into the second quadrant. So this boundary 
Umin = In tan w can be called an absorbing barrier^^. 
The vector passes positive w-axis counterclock-wise (see 
Sec. HIA). 

With all these, a mapping between our disordered QW 
model and the TBM with ODD is established regarding 
all the system parameters. Specifically, our Eqs. (32), 
(33) and (34) resemble Eqs. (21), (18) and (19) in Ref.^^, 
and the reflection and absorbing barriers are similar, too. 
Further borrowing the method in Sec. HI of Ref.^^, we 
directly find Ng 

Ng = -^-, With = 2 ((Intan^ ft)^) . (36) 

Using Eq. (32), we obtain the integrated DOS, 


iv,(w) = Ei-f 


4 In^ tan w 


and then the DOS, 


p{uj) = 


dNi 

dw 


1 

4 oj In^ Lo 


(37) 


(38) 


To conclude, we have shown that our disordered QW 
model possesses the physics of ODD. It is for this rea¬ 
son that, quite remarkably, the derivation of DOS for 
our QW model resembles to that in the original TBM 
with ODD2'^’2 i. To make this connection between our 
QW model and the TBM with ODD clear is the main 
contribution of this section. We highlight the two cru¬ 
cial steps: (i) linking the “spinor” components of the 
eigenstate through the transfer matrices, and (ii) the in¬ 
terpretation of the eigenstate as a vector moving in the 
complex plane when counting the number of states. 

The localization length for quasi-energies around 0 can 
be derived in a similar way^^ and the result is: 


cr^ In ui 

4 In^ tan w 4 In w 


(39) 


Equation (39) shows that the localization length diverges 
as OJ approaches 0, which is consistent with the previously 
mentioned fact that the state with w = 0 is delocalized. 
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D. Numerical analysis of the DOS 

The derivation of DOS in Sec. Ill C involves some ap¬ 
proximations, so we need numerical simulations to check 
the analytical results. Specifically, we use Eqs. (22), (29) 
and (30) to obtain the integrated DOS numerically, and 
then compare our numerics with the analytical expres¬ 
sion given by Eq. (37). Given one disorder realization 
and one quasi-energy w, we use the recursive relation in 
Eq. (22) to obtain and then it is substituted into 

Eq. (29) to obtain j, and finally we get Ni{uj) through 
Eq. (30). Note that a randomly chosen uj may not be 
an actual quasi-energy value associated with a particular 
disorder realization. However, if the system is sufficiently 
large, the quasi-energy values will cover the vicinity of 0 
quite densely. For this reason, a randomly chosen uj will 
not cause noticeable error in terms of the counting of 
states. 



FIG. 3. (Color online) Relation between integrated DOS and 
quasi-energy lo , shown via In (77/(w) — |) as a function of 
In I lntanta|. The QW chain is of size N = 3 x 10"*. The (red) 
solid line is from direct numerical calculations, the (blue) 
dashed line is a linear fit, and the (green) dash-doted line 
is our theoretical curve. The linear fit is applied to the do¬ 
main ln|lntanta| € [1,2], corresponding to the quasi-energy 
domain uj € [6.18 x 10“^, 6.60 x 10“^]. 

The analytical relation between Ni{uj) and w is given 
by Eq. (37). Alternatively, 

In ^A^/(w) — ~ —2 In I In tan w I. (40) 

Figure 3 depicts In (iV 7 (a;) — i) as a function of 

In I In tan w I to check this theoretical prediction. The 

2 

theoretical intersection on the y axis is In ^ « —1.40 
and the slope of the curve is —2. Our numerical results 
agree with theory well in the main domain of our inter¬ 
est. However, for w larger than « 0.066 (equiva¬ 
lently, In I In tan w I < 1), theoretical results deviate from 
the numerical data, implying the failure of the analyti¬ 
cal approximations made in Sec. HIC. This is expected 
as a too large oj leads to errors in Eq. (34) and then in 


Eq. (35). In the case of w < « 6.18 x lO”"* (equiv¬ 

alently, In I Intanwj > 2), the system size N is no longer 
large enough for a reliable statistical analysis, so the cor¬ 
responding numerical results also start to deviate from 
our theoretical predictions. 

E. A numerical study of the self-correlation of 
delocalized states 

Here we numerically check whether the average two- 
point correlation of a delocalized state with w = 0 de¬ 
cays polynomially. We use many realizations of disorder 
to obtain an average correlation function. This is differ¬ 
ent from our previous calculations where only a single 
realization of disorder is needed. Analytically, assuming 
that a dimensionless product of disorder strength and 
two-point separation is much larger than unity^^, the cor¬ 
relation exponent is shown to be —3/2. This theoretical 
prediction is checked here by use of Eqs. (17) and (18), 
which depicts the eigenstate structure of our disordered 
QW model. 



FIG. 4. (Color online) Dependence of correlation on the 
system size with the disorder strength fixed, as shown by 
In (|i/)(n)|^|i/>(l)|^) versus ln(n — 1), averaging over 2000 dis¬ 
order realizations. Here |i/>(n)p is the probability of the wave 
function at site n, and (|'i/)(n)p|'i/)(l)p) is the averaged two- 
point correlation, with one point fixed to be the site 1. From 
top to bottom, the system size is set to be A = 50, 100, 
400 and 2000 respectively, and the linear fitting curves have 
slopes —0.80, —0.98, —1.31 and —1.53. The disorder strength 
is fixed to be A = 0.4. 

In Fig. 4, the disorder strength is set to be A = 0.4, and 
the system size varies from iV-|-2 = 52 to 2002. When the 
two-point separation increases, the correlation exponent 
increases from 0.8 to 1.5 and stays almost stable at 1.5. 
Figure 5 shows how the correlation varies with the disor¬ 
der strength. The general observation is that increasing 
the disorder strength will increase the correlation expo¬ 
nent but the exponent again tends to saturate around 
—3/2. These numerical results are consistent with the 
early theoretical prediction of ODD^®’^^. However, we 






point out that if N and A are too large, the statistical 
fluctuations become more pronounced due to our limited 
number of realizations of disorder. 


o A = 0.2 
□ A — n /I 



denoted d. The initial state of the QW model is prepared 
with entries /3o = Q^i = l/VS and all other entries 0. It 
can be easily checked that this initial state is precisely 
the 0 -mode of the system (note that 6i is chosen to be 
tt/2). Then, we slowly reduce during the QW pro¬ 
cess, until 9n = Sn- To be more specific, the proposed 
adiabatic protocol can be achieved by introducing a slow 
time dependence to 9 in Eq. (7), i.e., 

9nit) = 9{t) + 5n, (41) 

with n G [1, N] denoting the bulk-site index, i5i = 0, and 
9{t) to be further specified below. 

The QW mapping operator U associated with 0„(t) 
is denoted as U{t). The initial state |'0(O)) is localized 
at the first two sites, with 1/(0) |i/;(0)) = |^/’(0)). The 
time-evolving state at time t is denoted \ip{t)), obtained 
by 


FIG. 5. (Color online) Dependence of correlation on dis¬ 
order strength with the system size fixed, as shown by 
In (|'i/)(n)|^|'i/'(l)P) versus ln(n— 1), averaging over 10000 dis¬ 
order realizations. \ip{n)\‘^ is the probability of wave function 
at site n. System size N + 2 = 202. Symbols circle, rectangle, 
and triangle represent A = 0.2, 0.4 and 1.0 respectively, and 
the linear fitting curves have slopes —0.86, —1.16, and —1.55. 


IV. EXPERIMENTAL PREPARATION OF THE 
O-MODE IN DISORDERED QW 

It is now clear that when the disordered local rotation 
angle variables fluctuate around zero (i.e., 0 = 0 in 
Eq. (7)), then the 0-mode (eigenstate with a; = 0) in 
our disordered QW model reflects the physics of ODD. 
However, if 0 7 ^ 0, then the corresponding 0-mode be¬ 
comes unrelated to ODD physics. For example, if 
slightly fluctuates around tt/2, then the 0 -mode will still 
be highly localized around the sites 0 and I, with negli¬ 
gible proportion in all other sites. 

The 0-mode with 0 = 0 is in general delocalized and 
hence it is hard to prepare in experiments. To address 
this issue, we note that the highly localized 0 -mode asso¬ 
ciated with 9 = IT 12 is a good starting point. We propose 
to connect this localized 0 -mode with our target 0 -mode 
possessing ODD physics by an adiabatic protocol^®”®^. 
That is, by slowly tuning the value of 9 from 7 r /2 to 0, 
we may reach our target 0 -mode from the localized 0 - 
mode. 

Consider then a conventional adiabatic evolution pro¬ 
tocol, through which the parameters in the QW opera¬ 
tor {7 are tuned slowly. Note, however, that the boundary 
rotation angles 9q and 6 *Ar+i must be fixed to ensure the 
conservation of probability inside the QW chain. An adi¬ 
abatic process reflecting this constraint is as follows. At 
first, the system is set as 9q = —tt/2, 9\ — 0jv_|_i = tt/2 
and 9n = 7 r /2 -|- with n G [ 2 , A] and being random 
angle fluctuations. The mean value of (5„ over A sites is 


m)) = U{t) • C/(t - 1) • • • U{1) ■ U{0) \m) ■ (42) 

For the sake of comparison between the time evolving 
state IV'(O) our target 0-mode state, we define the 
exact zero-quasienergy eigenstate of 17(f) as |V’°(^)) (with 
U{t) |'i/’°(f)) = |^/’°(f))). Numerically we can directly 

diagonalize U{t) to get |i/)°(f)). Our hope is to reach 
|'0°(f)) through the time evolving state \Tp{t)) emerging 
from our adiabatic protocol. Indeed, the adiabatic the¬ 
oremstates that \4’{t)) ~ if the adiabatic 

conditions are fulfilled. 

We have numerically simulated the process depicted 
in Eq. (42), and then compare with |'(/;°(f)). Their 

overlap probabilities | ('(/'(f)|V'°(f)) P versus t is plotted to 
check the performance of a certain specific protocol. In 
the following, by specifying 9(t) differently, we examine 
two protocols to realize the adiabatic process and hence 
the preparation of the target 0-mode state that reflects 
the physics of ODD. 

A. Tuning 6 at a constant rate 

In this case we decrease the bulk 9n at a constant rate 
with respect to the evolution time. Specifically, 9(t) in 
Eq. (41) is given by 

9(t) = 9(0) - rt, (43) 

where f = 0,1, 2,..., T is the evolution time, r = 9{0)/T 
is the constant decreasing rate, and 0(0) = tt/2. The 
obtained state fidelity | (^/>(f)|'0°(t)) p versus t is plotted 
in Fig. 6. 

Figure 6 shows that for some realizations of disorder, 
the fidelity near the final stage of the evolution decreases 
significantly. The difference seems to be related to 5, 
the actual mean value of the random fluctuations in 
a particular realization of disorder. In particular, the re¬ 
alization with 5 « —0.076 (green dash-dotted line) has 






9 


1 


0.9 


■+0 

^ 0.8 

¥ 0.7 


0.6 


0 100,. , 200 300 

time t 

FIG. 6. (Color online) Overlap probability between the 
actual time evolving state |V’(^)) s'-nd instantaneous 0-modes 
|V'°(t)) for 4 realizations of disorder in numerical experiments. 
The inset is a magnified view of the tail part. The (red) solid, 
(pink) dashed, (blue) dotted and (green) dash-dotted lines 
represent 4 different realizations of disorder with different 5 
(shown on the figure panel). The disordered chain has totally 
A'^-l-2 = 20 sites, with the disorder strength given by A = 0.7. 
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a final fidelity below 0.6. To understand this, we inves¬ 
tigate the gap between the 0 -mode and its neighboring 
mode, which is found to decrease with t. When 6{t) gets 
close to 0 , the 0 -mode is not well separated from the 
bulk modes, and the gap becomes quite small. Com¬ 
pared with other three realizations, the realization with 
S « —0.076 has a gap size of approximately half of oth¬ 
ers from t « 250 to 300, so this small gap has caused the 
most pronounced nonadiabatic transitions. To confirm 
this, we increase the total evolution time and indeed a 
better performance can be obtained (see Fig. 8 presented 
later). By contrast, for other realizations in Fig. 6 , the 
final fidelity is high (above 0.95), an indication of good 
performance due to the associated relatively large gaps. 
To summarize, the performance of this adiabatic proto¬ 
col is determined by the total evolution time T and the 
gap size in the final evolution stage. One can always im¬ 
prove the performance by increasing T. In contrast, the 
gap size is sensitive to the details of an actual realiza¬ 
tion of disorder. As an observation from our numerical 
results, cases with a negative <5 tend to have a smaller 
gap size around the final evolution stage than cases with 
a positive S. 


quasi-momentum. The gap between the bands is 29 at 
fc = 0. The 0-mode sits in the center of the band gap. 
We are thus motivated to design the following protocol 
by roughly assuming that the gap between the 0 -mode 
and the bulk spectrum is proportional to 9: 

^d(t) = -xm. (44) 

In this new protocol, the rate of change ^9{t) oc instan¬ 
taneous gap oc instantaneous 9{t) . As the gap decreases, 
the rate of change also decreases to keep the process be¬ 
ing sufficiently adiabatic. Therefore 9 is an exponential 
function of t, 

9{t) = 9{0)e-^\ (45) 

where A is the exponential decay rate of 9. Using this 
protocol, 9n can be explicitly expressed as a function of 
t: 


n = 0 , 

n = 1 , , , 

^d(T) + 5„ ne[2,iV], ^ ^ 

n = N + 1. 

Here j^^9(T) is to make sure that ~ 

J2n=i at the final time t = T. Note also that at site 
n = 1, 6*1 (0) = 6*(0) = 7 r/ 2 , which ensures that the initial 
0-mode is the exact eigenstate of the QW propagator at 
time zero. 

Figure 7 shows the performance of this protocol. For 
positive 6, the overlap probability at final time is quite 
high (above 0.998). Interestingly, similar to the previ¬ 
ous protocol in which we sweep 6 * at a constant rate, 
the fidelity degrades in cases of J < 0. In addition, in 
some realizations of disorder, the gap size may be erratic 
during the last stage of the adiabatic protocol, especially 
when S turns from positive to negative. This explains the 
relatively poor performance for the case with S = —0.108 
in Fig. 7. 

Nevertheless, we can further improve the fidelity by 
increasing the total evolution time T or decreasing A in 
our exponential protocol. Panel (a) of Fig. 8 how fidelity 
changes with T. As a comparison, in panel (b) of Fig. 8 we 
show the parallel fidelity vs T if d is swept at a constant 
rate. It is seen that overall, tuning 9 exponentially as 
is done here is much better than tuning 9 at a constant 
rate. 



B. Tuning d exponentially 

To understand our motivation of this alternative pro¬ 
tocol, we first discuss the gap size of the clean system, 
where the bulk is uniform (i.e., = ^)- In this case, 

two quasi-energy bands emerge and the dispersion re¬ 
lation is given by cosw = cos 6 * cos /c’'*, where k is the 


C. Correlation exponents in numerical experiments 

We have shown in the previous subsection how to pre¬ 
pare the 0-mode state possessing the physics of ODD. 
Here we aim to show that states prepared in this manner 
can indeed manifest the correlation exponent character¬ 
istic of ODD physics. In doing so we need to perform 
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FIG. 7. (Color online) Overlap probability between the 
actual time evolving state IV'(^)) a-nd instantaneous 0-modes 
|i/’°(l)) for 2 different types of disorder realization. The chain 
has -b 2 = 20 sites, total evolution time T = 90, disorder 
strength A = 0.7, and the parameter in the exponential pro¬ 
tocol is characterized by A = 0.0562. (red) Circles are for 
a case with the averaged angular disorder 5 = 0.064 being 
positive, with the overlap probability above 0.998 at the final 
time. The inset shows more details, (blue) Triangles for a 
case with the averaged angle disorder 5 = —0.108 being nega¬ 
tive. In this case, the final overlap probability is only around 
0.65, which means that this protocol is still not working well 
with T = 90. 


averaging over many realizations of disorder. We use the 
exponential adiabatic protocol in our numerical experi¬ 
ment. To benchmark our numerical experiments, we also 
analyze the correlation exponent using the exact delocal¬ 
ized 0-mode state obtained from Eqs. (17) and (18). 

Before presenting our results, we first discuss two mi¬ 
nor issues. The first is related to the fact that the spinors 
represented in Fig. 1 involve two different sites. That is. 
In a real experiment, what is measured is likely the prob¬ 
ability at each site, whereas in our analytical study, we 
treat (/3n-i oin)^ as one “spinor”. However, we find that 
this difference has little effect on the correlation expo¬ 
nent. The other issue is that we have fixed 6*i to be 7r/2 
(hence not random) (see Sec. IV for details). Again, it is 
checked that this does not affect our analysis. 

We also note that the —3/2 correlation exponent was 
derived under the assumption that the product of the di¬ 
mensionless disorder strength and two-point separation 
is much larger than unity^^. In real experiments, the 
QW chain might not be long, so we are limited to rela¬ 
tively small two-point separation. That means we should 
choose strong disorder strength to fulfill this assumption. 
Figure 9 presents our results from numerical experiments 
based on an exponential adiabatic protocol starting from 
a highly localized state, as compared with a direct inves¬ 
tigation using the exact delocalized 0-mode states. For 
two different chain length, the two-point correlation ex¬ 
ponents in our numerical experiments are found to be 




FIG. 8. (Color online) Overlap probability versus t for dif¬ 
ferent protocol duration T, for an exponential protocol (a) 
(Eq. 45)) and the previous constant-rate protocol (6). In both 
protocols, the disorder realization is the same as the one with 
S = —0.1083 in Fig. 7, and A-l-2 = 20, A = 0.7. (a) From top 
to bottom, T equals 240, 210, 180, 150, 120 and 90. The corre¬ 
sponding values of A is chosen to be A = — ln(0.01/(7r/2))/T. 
(b) From top to bottom, T equals 400, 300, 240, 180 and 
90. In both panels, a larger T results in a better fidelity of 
the final state. However, the exponential protocol in general 
requires less time to achieve the same fidelity. 


— 1.48 and —1.36, as compared with —1.6 and —1.5 ob¬ 
tained from pure theory. Certainly, the agreement be¬ 
tween these two sets of data can be further improved if 
we further increase T. The conclusion is that our adi¬ 
abatic protocol applied to our disordered QW model is 
also useful in the actual demonstration of the two-point 
correlation characteristic of ODD physics. 

For small systems with weak disorder, the analytical 
correlation exponents are not available^^. To motivate 
experimental studies on this matter, below we further 
exploit our setup to investigate how the two-point corre¬ 
lation changes with weak disorder strength A and system 
size (N+ 2). 

We choose 4 different system sizes with a fixed and 
weak disorder strength A = 0.4. In particular, we let 
N + 2 = 12, 22, 32 and 42. The results are shown in 
Fig. 10. For each case, we show statistical results ob¬ 
tained from analytical treatment of the 0-mode with dis¬ 
order and from our exponential adiabatic protocol that 
starts from an initial localized state. The results obtained 
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FIG. 9. (Color online) Correlation function 

In n)p|i/)(r, 1)1^) versus ln(n — 1), averaged over 1000 

disorder realizations, (a) disorder strength A = 1, system 
size + 2 = 32; (6) A = 1, A + 2 = 42. The total evolution 
time T is chosen to assure satisfactory fidelity in the adia¬ 
batic preparation of the 0-mode, with T = 400 in panel (a) 
and T = 600 in panel (6). In both panels (red) circles denote 
results from solving the 0-mode analytically; whereas (blue) 
stars denote results obtained from our adiabatic preparation 
of the 0-mode with the exponential protocol. Solid line and 
dash-dotted line are the associated linear fitting curves over 
a regime without much fluctuation. The slopes of the fitting 
curves reflect the correlation exponents. 

from such two totally different methods agree very well 
because they yield almost the same slopes from the fit¬ 
ting straight lines, for all the four cases shown. The good 
fitting by the straight lines indicates a polynomial behav¬ 
ior of the two-point correlation function, but now with 
correlation exponents given by —0.447, —0.588, —0.645 
and —0.769, for N = 10, 20, 30 and 40, respectively. 
These exponents are far from -3/2, but shows a tendency 
to approach -3/2 as the system size increases. Further 
increasing the value of A also increases the magnitude of 
the correlation exponent. These results should be of ex¬ 
perimental interest as well and invite further theoretical 
developments in studies of the physics of ODD. 

V. SUMMARY 

To summarize, we have shown that the physics of ODD 
can be investigated by a disordered QW model. The as- 



FIG. 10. (Color online) Correlation functions with weak 
disorder in a QW model, shown via In (|'i/)(T, n)|^|i/)(r, l)p) 
versus ln(n— 1), averaged over 1000 disorder realizations. The 
total evolution time T is chosen to make sure the adiabatic 
protocol can yield a satisfactory fidelity of the 0-mode state. 
For example, if A or A is increased, T is increased also (See 
Sec. IV B). Here A = 0.4, and from top to bottom, the system 
size is A -I- 2 = 12, 22, 32 and 42 respectively. The slopes of 
the curves fitting the results using the exact 0-mode (red solid 
line) are —0.45, —0.59, —0.65 and —0.77, whereas the slopes 
of the curves fitting the results arising from our adiabatic 
protocol (blue dash-dotted line) are —0.42, —0.59, —0.66 and 
—0.73 respectively. The symbols and the lines share the same 
meaning with those in Fig. 9. 

sociated exotic features in the delocalization and in the 
wavefunction correlation are derived and numerically ver¬ 
ified. Because the physics of ODD is rarely cleanly ob¬ 
served in actual experiments, our results will possibly 
motivate ongoing QW experiments as a new platform to 
study the physics of ODD. To facilitate such efforts, we 
proposed and analyzed adiabatic protocols to prepare the 
exotic delocalized 0-mode state with good fidelity. Our 
numerical experiments show that the delocalized 0-mode 
states thus obtained can directly show the correlation ex¬ 
ponent -3/2 in the regime predicted by existing theory. 
Our numerical experiments also show that much differ¬ 
ent correlation exponents emerge if the product of the 
system size and the disorder strength is relatively small. 

Appendix A: From Eq. (15) to Eq. (19) 

Here we show how to derive Eq. (19) from Eq. (15). 
Multiply both sides of Eq. (15) with , and decom¬ 
pose Tn using the following identity 

^ fe^‘^sec8n - tan6i„ \ 

" \^ —tan0„ e~^‘^ sec 9nJ 

_ ® ^ sec9n — tan0„\ 0 \ 

“ 0 e“®T J i^-tan0„ sec0„ / V ^ / ’ 

(Al) 
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then Eq. (15) becomes 



/e*2 0 \ / sec0jv — tan0Ar\ 

0 J \^—tan^AT sec0jv ) 

fe''^ ® ^ sec0jv_i —tan^AT-iA 

0 e“*‘^y \^—tan6>7v_i sec0Ar_i y 

/ sec 01 — tan0i\ ® ^ ^ ^ 

\^—tan 01 sec 01 / \ 0 e“®T y ye“*Ty ’ 

(A2) 


Replace 

titles 




in Eq. (A2) with the iden- 



Eq. (A2) then becomes 



/ sec0„ — tan0„\ 0 \ 

\^—tan0„ sec0„ / V 0 e“®“y 




Multiply matrix P~^ from the left of both sides of 
Eq. (A4) and insert the identity / = P~^P between 
neighboring matrices in the right hand side, where P~^ = 
P = {(^x + we will arrive at Eq. (19) because 


p_i 0 ^ p _ cosw isinw\ 

0 e~^^j ~ l^isina; cosw j ’ 


p_x f sec On 

I — tan0„ 


— tan0„\ _ Aani9„ 0 \ 

secOn ) ~ 0 tani?„y ’ 


(A5) 


where = f - 


Appendix B: More on the boundary conditions 

Previously we employ one specific boundary condi¬ 
tion to study the physics of ODD, but leave three 
other boundary conditions unexplored. Here we will 
briefly summarize the special quasi-energies and the cor¬ 
responding states'^®’®' for these different boundary con¬ 
ditions. Given the bulk 0„ = 7r/4 -|- (5„ with |(5„| < 7r/4, 
then the boundary condition (0o, On+i) = (-7r/2,7r/2) 
[(7r/2, —7r/2)] will lead to the edge states with quasi¬ 
energy Lo = 0 or TT localized around the boundary site 
n = 0 [n = TV -I- 1]. For convenience, we assume = 0 
in our qualitative discussions below. 

Interestingly, the 0 or tt quasi-energy states are absent 
under the boundary conditions (0o, 0Ar+i) = (7r/2,7r/2). 
For the case of (0o,0Ar+i) = (—7r/2, —7r/2), it can be 


shown that there exist localized edge states with quasi¬ 
energies slightly differing from 0 or tt. These features 
are also relevant to understand the topological properties 
in Here we elaborate these features using the 

transfer matrix formalism (TMF). Following the same 
method in Sec. HI, the relation between 2 boundaries 
given by Eq. (15) can be written in the form analogous 
to Eq. (19): 


-P -(“), (Bl) 

\^oy i^zsinw cosw y V*/ ^ ^ 

(“)=cy™“ ®“A.p.(!). (B2) 

yt J yzsinoj coscu J yOJ ^ ^ 

Here = f — ^ and P is given in Eq. (19). Eq. (Bl) is 
for the boundary condition (0o,0Ar+i) = (7r/2,7r/2) and 
Eq. (B2) is for (0o,0Ar+i) = (-7r/2,-7r/2). 

In the case of Eq. (Bl) and using the same language as 
in Sec. HI A, an actual quasi-energy w needs to bring a 

vector initially at the y-axis, to the x-axis, . For 

simplicity, we assume the vector goes from the positive y- 
axis to the negative x-axis. w = 0 or tt certainly cannot 
accomplish this task since it will let the vector stay in 
y-axis. Let us check if a small value e which slightly 
above 0 can be the quasi-energy, using Eq. (22) with 
On = 7r/4, (/)i = -K12, (pN = TT — e and = tt/S. It then 
follows that tan cpn should approach 0 from —c» (that is, 
after the vector enters the second quadrant). However, 
this cannot be true since cot^(7r/8) >> 1 will prevent 
tan (pn from approaching 0. Together with other simple 
considerations, it is seen that under the above boundary 
condition, a; = 0, tt and any value near them cannot be 
the quasi-energies of the system. 

In the case of Eq. (B2) , the vector should go from the 
x-axis to the y-axis. For simplicity, we assume the vector 
goes from the positive x-axis to the positive y-axis. This 
corresponds to tan (pn going from 0 to c». It is obvious 
that w = 0 or TT cannot achieve this goal. Again we con¬ 
sider a small value w = e. Now the factor cot^(7r/8) >> 1 
in Eq. (22) will speed up this process, thus indicating that 
a small oj = e may satisfy Eq. (B2). In addition, accord¬ 
ing to Fig. 2, when cp is smaller than 7r/4, the length of the 
vector tends to decrease exponentially, and after it passes 
7r/4, the length starts to increase exponentially. There¬ 
fore, the corresponding eigenstate is sharply localized at 
both edges. Except for this particular e, we may expect 
that a vector with a slightly larger oj may pass two more 
quadrants to reach the negative y-axis such that it can be 
another quasi-energy of the system. But this is not true 
because the vector cannot goes from the positive y-axis 
to the negative x-axis. Hence, this small quasi-energy e is 
well-separated from other quasi-energies. Until a quasi¬ 
energy uj becomes large enough to cross the 2nd quadrant 
(i.e., from the positive y-axis to the negative x-axis), no 
other Lu can satisfy Eq. (B2). 
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Appendix C: Other special quasi-energies in the 
disordered QW 


Obuse et al numerically showed that w = ± 7 r /2 can 
be also special quasi-energy values with singular DOS, 
which hence indicate the presence of ODD in disordered 
QW. Here we use the method developed in Sec. Ill to 
discuss these special quasi-energy values.x 

We start with Eqs. (11) and (14) in Sec. II A. Without 
loss of generality, we choose uj = 7 r/ 2 . Then the chain 
relation analogous to Eq. (15) will be 


Cn 


— sin 9o 


N 


n -^0 (sin 0^+1 ' 

n—1 ^ ' 


(Cl) 


= ia, sec tan 9„. 


Define 


In our set-up, TV -|- 2 is the total number of sites in the 
disordered QW chain (See Fig. 1). Each bulk site cor¬ 
responds to one transfer matrix, and totally N transfer 
matrices are involved in the calculation. When N is odd. 


Boundary condition 

u! = ±^, N even 

ta = ±y, A odd 

00 = ^ = 0iV+l 

Y 

N 

n TT /n TT 

^0 — — J, ^N + 1 — J 

N 

Y 

n TT /n TT 

(70 — f7iV-|-l — —-2 

N 

Y 

o 

II 

1 

+ 

Y 

N 


TABLE I. The existence (Y) or nonexistence (N) of modes 
under different boundary conditions. In the bulk, values of Qn 
(1 < n < A) are assumed not to satisfy 7r/4 — 0^/2 = j • 7r/2 
(j is an integer). 


one transfer matrix will be left if we pair those transfer 
matrices according to Eq. (C2). This leads to 


Pm = T2m • 22m-1, 


(C2) 


so 


Pm = (tan 6 » 2 mtan 6 ' 2 m-l - sec 6 » 2 m sec 6 > 2 m-l) - 2 + 
(sec 02 m tan 02 m-l - tan 02 m sec 02 m-l) • CTj,. 

(C3) 

Expressing Pm in the basis of cxy, we have 
'-cotrT2mtanrT2m-i 0 


p — 
J m, — 


0 


- tan 92m cot ■ 02 m-l J ’ 
(C4) 


where = f — So in the ay basis for even TV, 


N 


Hr,. 

n—1 


x+ 0 

0 A_ 


(C5) 


= 2 ~’ rrs'i 

n—1 ^ ' 

[(A+ -I- A'_) • / -I- (A+ - A'_) • ay] , 

where A(^ and A(_ are obtained from Eq. (C6) by substi¬ 
tuting N with N—1. Different from the case of even 
TV, the additional ax and tr^ flip the eigen spinors of 
CTy, resulting in the opposite conclusions. In particu¬ 
lar, boundary conditions 9 q = 9 m+i = — 7r/2 or 

00 = —'K12, 9^+1 = 7r/2 will give rise to a; = 7r/2, while 
00 = Qn+i = ±7r/2 cannot. 

We summarize the results in the Tab. I. Those states 
with exactly quasi-energy ±7r/2 are delocalized. For ex¬ 
ample, in the case of even TV and 0o = 0jv+i = ~7’‘/2, we 
substitute Eq. (C7) into Eq. (Cl) and get 


CAf 1 „■ I = *CoA+ 


(C9) 


with 

1 TV 

A+ = Al = (—1)"^ cot'0Artan'0jv-i • • • cot'02 tan-di. 

(C6) 

Returning to the basis, we have 

^ 1 

n = - [(A+ + A_) • / + (A+ - A_) • ay] . (C7) 

n—1 

We substitute Eq. (C7) into Eq. (Cl) and find that 
the boundary conditions 0o = 0Ar+i = ±7r/2 will make 
Eq. (Cl) hold, while 0o = 7r/2, 9n+i = —7r/2 or 
00 = —'K12, 0Ar+i = 7r/2 cannot. This conclusion is in¬ 
dependent of the actual values of 0„ (n = 1,2 • • • TV), so 
whether w = 7r/2 is the quasi-energy of the system is 
determined by the boundary conditions, as well as the 
parity of the number of system sites. 


Therefore, the spinors at two boundaries are the eigen 
spinor of ay, and they are connected by A+ in Eq. (C 6 ). 
In general A+ « I because cot and tan-dfc (j, k S [I, N] 
are arbitrary indices) will approximately cancel each 
other given that 0 j\j, are drawn randomly from a given 
distribution. This resembles the 0-mode in Sec. IIB. Note 
that, the delocalized 0 -mode requires 0„ to be drawn 
from a distribution symmetric with respect to 0 = 0 (we 
choose 9n G [—A, A] in our study), whereas the delocal¬ 
ized ± 7 r /2 states do not have this constraint. However, 
the advantage of a delocalized state at w = 0 is that it can 
be obtained from localized w = 0 state through an adia¬ 
batic protocol (See Sec. IV). By contrast, the uj = ± 7 r /2 
states cannot be obtained in this manner. The reason is 
simple. States with oj = ± 7 r /2 are delocalized regardless 
of 0 , the mean value of 0 „; whereas a delocalized ui = 0 
state requires 0 « 0 . 
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